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Abstract 

Depending on temperature the modified Hodgkin-Huxley (MHH) equations exhibit a va- 
riety of dynamical behavior including intrinsic chaotic firing. We analyze synchronization in 
a large ensemble of MHH neurons that are interconnected with gap junctions. By evaluating 
tangential Lyapunov exponents we clarify whether synchronous state of neurons is chaotic or 
periodic. Then, we evaluate transversal Lyapunov exponents to elucidate if this synchronous 
state is stable against infinitesimal perturbations. Our analysis elucidates that with weak gap 
junctions, stability of synchronization of MHH neurons shows rather complicated change 
with temperature. We, however, find that with strong gap junctions, synchronous state is sta- 
ble over the wide range of temperature irrespective of whether synchronous state is chaotic or 
periodic. It turns out that strong gap junctions realize the robust synchronization mechanism, 
which well explains synchronization in interneurons in the real nervous system. 
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■ In the rat hippocampus interneurons show the high frequency synchronization during the gamma 

oscillation (^40Hz) and sharp wave burst (^200Hz) [1], and such simple synchronization of a 
lO ' large ensemble of neurons has attracted much attention of theoretical researchers [2-10]. One 

major analysis for these studies is the phase reduction method [9-11], in which phase vari- 
ables are utilized to represent the periodic behavior of neurons. The phase reduction method is, 
however, applicable only to the case of infinitesimal interactions. Moreover, if neurons behave 
aperiodic, phase variables are indefinable. The general synchronization properties of strongly 
coupled neurons thus remained unclear, especially in the case of chaotic neurons. 
• ■ Meanwhile, studies of synchronization of a large ensemble of chaotic oscillators have made 

r> \ a remarkable progress in recent years [12-15]. The major targets of these studies are simple 

chaotic oscillators such as Lorenz equations and logistic maps. Synchronous state of these os- 
cillators is characterized by two types of Lyapunov exponents: tangential Lyapunov exponents 
and transversal Lyapunov exponents. While tangential Lyapunov exponents clarify whether syn- 
chronous state is chaotic or periodic, transversal Lyapunov exponents elucidate if synchronous 
state is stable against infinitesimal perturbations. In the present study, we employ these sophisti- 
cated techniques in chaos synchronization theory to investigate synchronization of neurons. We 
show that tangential and transversal Lyapunov exponents enable us to analyze stability of syn- 
chronization in a large ensemble of neurons for arbitrary neuron dynamics and arbitrary strength 
of interactions. 

The concrete target of the present analysis is a network of iV(> 2) spiking neurons that 
obey the modified Hodgkin-Huxley (MHH) equations [16, 17]. The MHH equations are four- 
dimensional nonlinear differential equations, which include temperature-dependent scaling fac- 
tors p = a\ T ~ Tq)/W and = " T ° )IW (See Ref. [17].) With T changing, a MHH neuron shows 
a variety of dynamical behavior including chaotic firing as shown in Fig. ^a) (ISI and so on 
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will be explained later.) For the sake of simplicity, we denote this MHH neuron dynamics by 
dx/dt — F(x) with neuron state vector x = (v, w\, . . . ,w„_i) T , where v represents the membrane 
potential and {w/} describe gating of ion channels. We assume that N neurons {x,} are intercon- 
nected with all-to-all gap junctions. Since gap junctions induce electric currents proportional to 
potential difference between neurons, the dynamics of the neural networks is expressed as 

dx- 

— i=F(x / ) + (/ / /c,0,...,0) T , i=l,...,N (1) 
at 

with 

li = Jf E ( V J ~ v = # £ (*fl - *n ) ' ( 2 ) 

where constant c = 1.0 /jF/cm 2 is capacitance of the membrane and /, is the electric current 
induced by gap junctions. 

The above dynamics can be generalized to the form 

^= F ( x + flG( Xi , X ,-). (3) 

Therefore, we investigate synchronous state in this general mean-field dynamics. We assume 
stationary synchronous state x\ = . . . = x* N = x* , which obeys 

dx* 

= F(x*)+gG(x*X). (4) 

dt 

To elucidate the stability of this synchronous state we investigate perturbed state x, = x* + Sx,-. 
We define Jacobi matrices such that F(x* + Sx) = F(x*) + F'(x*)8x + (higher order) and G(x* + 
5x i , x* + 8x 2 ) = G (x* , x* ) + G J (x* , x* )Sxi + G' 2 (x* , x* ) 8x 2 + (higher order) . Then, Taylor series 
expansion to the first order yields 

= [F'(x*)+gGUx*,x*)]8x, 

+gG 2 (x*,x*)l£Sx 7 , (5) 

The naive evaluation of this A^-body dynamics brings about an eigenvalue problem with the 
large size of matrix. We hence define mean state x = (1 /Af)£ ; x,- and obtain the dynamics of its 
deviation Sx = (l/N) 8x, in the closed form 

For this n-dimensional linear dynamics, we can define the spectrum of n Lyapunov exponents {A,J}, =1 
These exponents are the so-called tangential Lyapunov exponents. To the first order, Eqs. © and 
(|6} are equivalent to 

— (x* + Sx) = F(x* + Sx) + 8 G(x* + 8x,x* + Sx). (7) 
dt 

Solving Eqs. (0J and Q numerically we can calculate time evolution of sufficiently small de- 
viation Sx. Evaluating this time evolution of Sx by the well-known computational method for 
Lyapunov exponents [18] we can calculate {A.J} numerically. Note that in this calculation of {A-}} 
we do not have to solve the huge A^-body dynamics in Eq. (|3j- Since replacement of x* in Eq. (|4} 
by x* + Sx gives the same dynamics as Eq. 0, {X\} indicate the characteristics of synchronous 
state x*, that is, when synchronous state x* is periodic (chaotic), the largest tangential Lyapunov 
exponent X\ takes the zero value (a positive value). 

We have evaluated mean state x by tangential Lyapunov exponents {kj}- We now investigate 
deviations around mean state: x, = x + Sx,-. Subtracting Eq. from Eq. l|3 we obtain the 
dynamics of Sx, in the closed form 

^=[F'(x*)+gGUx*,x*)]Sx,-. (8) 
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When synchronization is stable, all deviations {8x,} must converge to {0}. Since these N dy- 
namics of deviations are completely identical, it suffices to evaluate only one dynamics among 
them. For «-dimensional linear dynamics in Eq. (jSJi we can define the spectrum of n Lyapunov 
exponents {Xj} l=l ■ These exponents are the so-called transversal Lyapunov exponents. The 
largest transversal Lyapunov exponent Xj takes a negative value when synchronous state is stable 
in the sense of Milnor [14, 19]. To the first order, Eqs. and (jSJi are equivalent to 

j ( {x* + Mi) = F(x* + Si,-) +gG{x* + Sxi,x*). (9) 

Applying the computational method for Lyapunov exponents [18] to Eqs. and l|9), we can 
calculate {Xj} numerically. 

Let us apply the above analysis to investigating synchronization in networks of MHH neurons 
defined by Eqs. Q and ©. First, we calculate synchronous state x* in Eq. @. Note that in the 
present system the interaction term gG(x*,x*) in Eq. @ vanishes because of Eq. 0- Therefore, 
the behavior of x* in Eq. l|4} is completely the same as that of an isolated single MHH neuron. 
For the rough illustration of a single MHH neuron behavior, we define the k-th spike timing t (k) 
by the time when membrane potential v* = xj crosses the threshold value = —20 mV from 
below, and then calculate interspike intervals (ISIs) t(k+ 1) — t (k) (k — 1,2, . . .) in Fig. [TJ a )- 
Below T = 6.8 °C, ISIs take the single value around 650 msec, implying the simple periodic 
firing in which only one spike arises during the period. At T = 6.8 °C, however, period doubling 
bifurcation occurs so that the neuron fires twice during the period. After that, following typical 
period doubling cascade, the MHH neuron dynamics reaches the chaotic regime beyond T = 
7.3 °C, where ISI distribution becomes blurred. In this chaotic regime, we, however, observe 
several periodic windows, in which the behavior of neuron becomes periodic abruptly. 

Second, from Eqs. © and @, we calculate the tangential Lyapunov exponents {A-J} for the 
exact characterization of synchronous state x* illustrated in Fig.Qa). In the present system, {X\} 
are independent of g because of the interaction in Eq. (|2j. In Fig.^b), we plot the largest tan- 
gential Lyapunov exponents X\ as a function of temperature T. When synchronous state x* in 
Fig.|H a ) is periodic, X\ takes the zero value. In chaotic regime beyond T = 7.3 °C, X\ takes a pos- 
itive value, though we observe the several valleys of X\ corresponding to the periodic windows 
observed in Fig.Qa). 

Third, we calculate the transversal Lyapunov exponents {Xj} from Eqs. @ and {9j. {Xj} de- 
pend on parameter g. In Fig.^c) we calculate the largest transversal Lyapunov exponent Xj as a 
function of temperature T forg = 0.02 mS/cm 2 . When Xj takes a negative value, synchronization 
of MHH neurons can occur. 

When we assume weak gap junctions as in Fig.Qg = 0.02 mS/cm 2 ), the condition for syn- 
chronization of MHH neurons is rather complicated. Periodic synchronous state is stable in some 
conditions and unstable in other conditions. We also see stable chaotic synchronous state in some 
values of temperature T. Around T ~ 12 °C we find unstable periodic synchronous state inside 
the periodic window (X\ = and < Xj at T = 1 1 .9 °C) while we find stable chaotic synchronous 
state outside the periodic window (0 < X\ and Xj < at T = 12.1 °C). Actually, the numerical 
simulations of 100 MHH neurons in Fig|2]show the good agreement with the results of our analy- 
sis. Around T ~ 9.5 °C, however, synchronous state is stable both inside and outside the periodic 
window. With weak gap junctions the condition for synchronization is so complicated that its 
intuitive explanation is difficult. 

On the other hand, with strong gap junctions, synchronization of MHH neurons is stable 
over the wide range of temperature T . In Fig. [3] we plot the largest transversal Lyapunov ex- 
ponent Xj as a function of g for various values of temperature T. Eqs. and (jSJl show that 
when parameter g takes the zero value, transversal Lyapunov exponents {X^} take the same value 
as tangential Lyapunov exponents {Xj}. Therefore, synchronous state in Fig.[5]is periodic for 
T = 6.5, 7, and 11.9 ( °C) and chaotic for T = 7.5, 11, and 12.1 ( °C). In these six temperature, 
the behavior of MHH neurons are quite different from one another. However, all Xj take negative 
values if we increase the strength of gap junctions beyond g = 0.05. In all the temperature we 
investigate (5 °C < T < 15 °C), we find a certain value of g beyond which Xj always take a nega- 
tive value. Irrespective of whether synchronous state is chaotic or periodic, strong gap junctions 
induce synchronization of neurons. 
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In summary, we have studied synchronous state of a large ensemble of modified Hodgkin- 
Huxley (MHH) neurons assuming gap junctions among neurons. For the general mean-field 
dynamics in Eq. Q, we have evaluated N x n-dimensional deviation to define tangential Lya- 
punov exponents {A,^} and transversal Lyapunov exponents {Xj}. In Fig.^b), we have inves- 
tigated characteristic of synchronous state of MHH neurons by the largest tangential Lyapunov 
exponent X\. In Fig.Qc), we have elucidated stability of this synchronization by the largest 
transversal Lyapunov exponent X\. With weak gap junctions g = 0.02 mS/cm 2 , stability of syn- 
chronization of MHH neurons shows rather complicated change with temperature as shown in 
Fig.^c). However, in Fig. [3] and so on, we have found that with strong gap junctions, syn- 
chronous state is stable over the wide range of temperature. The strong gap junctions induce 
synchronization both in periodic and chaotic neurons, and that implies a pivotal role of gap junc- 
tions in synchronization of the large number of neurons 

It should be emphasized that the computational cost for Lyapunov exponents of four-dimensional 
MHH equations is not much higher than that of the three-dimensional Lorenz equations. Even 
when we assume dozens of ion channels in neuron dynamics, we would be able to calculate 
both Lyapunov exponents within the acceptable computation time. When synchronous state is 
periodic and interactions are infinitesimal (g <C 1), one can use the phase reduction method. 
In the present system of MHH neurons, however, A,j in Fig. [5] fluctuates dramatically as g in- 
creases, even if synchronous state is periodic. Since the synchronization property with finite g 
is quite different from that with infinitesimal g, calculation of transversal Lyapunov exponents 
is indispensable in investigating the present system. In the real nervous system many types of 
neurons are interconnected with chemical synapses. Some authors model the dynamics of chem- 
ical synapse in the manner as /, = {g/N)Y.j{vrev —Vi)sj, where constant v rev denotes reversal 
potential and sj obeys the dynamics dsj/dt = aF(v;)(l — s,-) — fisj with the sigmoidal function 
F(v) = 1/(1 + exp[— (v — 9. svn )/2]) [2,5,22]. In this case, synchronous state x* depends on g 
since G(x*,x*) does not vanish. Moreover, Jacobi matrices of G(x,-,Xj) are not constant but 
depend on x, and x ; -. Our analysis is applicable also to such complicated neural networks since 
their dynamics are written in the form of Eq. 0. Although pulse-coupled neural networks based 
on threshold-crossing spike timing cannot be written in the form of Eq. Q, we can employ the 
similar analysis by carrying out the decomposition of linear stability discussed in our previous 
study [8]. In that study, we have evaluated two types of Floquet matrices to show that periodic 
synchronous state of integrated-and-fire (IF) neurons are stable with only inhibitory chemical 
synapses. Interestingly, in the real nervous system, interneurons are found to be connected with 
inhibitory chemical synapses and gap junctions [20]. It turns out that networks of interneurons 
take the extremely ideal structure to induce synchronization of an ensemble of neurons. The 
present approach of stability analysis is applicable to a wide class of stability problems in neu- 
ral networks. Retrieval state in associative memory neural networks of spiking neurons can be 
investigated by the similar stability analysis [21]. More complicated neural networks including 
pyramidal neurons and interneurons [23] would also be analyzed by the present approach. 
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Figure 1 : (a) The interspike intervals (ISIs) in stationary state of an isolated single MHH neuron 
are plotted as a function of temperature T. The ISI distribution of synchronous state of N neu- 
rons x* is the same as a isolated single neuron since the interaction term in Eq. @ vanishes in 
the present system, (b) The largest tangential Lyapunov exponent X\, which characterizes syn- 
chronous state x* described in (a), is plotted. In the present system, X\ is independent of g. (c) 



The largest transversal Lyapunov exponent X\ is plotted for g = 0.02 mS/cm 
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Figure 2: The results of numerical simulations of 100 MHH neurons with g = 0.02 mS/cm 2 are 
plotted for (a) T = 11.9 °C and (b) T = 12.1 °C. Each dot represents spike timing in stationary 
state realized after t = 1.0 x 10 7 msec. 
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Figure 3: The largest transversal Lyapunov exponent A,j is plotted as a function of g for various 
values of temperature T. The numbers in the figure indicate temperature T. 
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